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^ r-| Abstract 

> : 

^ . We compute the nucleosynthesis bounds on the masses of stable Dirac and Majorana neu- 
trinos by solving an evolution equation network comprising of all neutrino species which 
in the Dirac case includes different helicity states as separate species. We will not commit 
ourselves to any particular value of the nucleosynthesis bound on effective number of light 
neutrino degrees of freedom N u , but present all our mass bounds as functions of AN U . For 
example, we find that the excluded region in the mass of a Majorana fi- or r- neutrino, 
0.31 MeV < mff < 52 MeV corresponding to a bound AN U < 0.3 gets relaxed to 0.93 MeV 
< mf, 1 < 31 MeV if AN U < 1.0 is used instead. For the Dirac neutrinos this latter constraint 
gives the upper limits (for T QCD = 100 MeV): m Vfi < 0.31 MeV and m Vr < 0.37 MeV. Also, 
the constraint AN U < 1 allows a stable Dirac neutrino with mP > 22 MeV. 



1 On leave of absence from the Department of High Energy Physics (SEFL), University of Helsinki, Finland. 



1 Introduction 



Primordial nucleosynthesis considerations have become a widely used tool to obtain limits 
on particle properties such as masses, couplings and lifetimes. Nucleosynthesis bounds arise 
from the tight agreement between primordial abundances of the light elements deduced from 
the observations and the theoretical predictions based on the standard big bang nucleosyn- 
thesis model (SBBN) Typically, any extension of the standard model, such as admitting 
large neutrino masses, could destroy the agreement between the theoretical prediction and 
the observational evidence. 

One of the quantities already studied in the context of nucleosynthesis is the mass of 
a stable neutrino Nucleosynthesis is sensitive to neutrino masses in the interval 

m u ~ 0.1 — 50 MeV, and the actual values of the bounds depend on the particular value 
adopted for the nucleosynthesis bound on the effective number of neutrino degrees of freedom 
AN U . The nucleosynthesis bound on AN U has proven to be difficult to pin down accurately 
and has been under constant revision over the past years Jl|, [j5]-|pl|. Most recently, doubts 
have been raised regarding the consistency of the standard big bang nucleosynthesis model 
with 3 massless neutrinos [TT], inducing a closer look into the issue of possible systematic 
errors Hl2| in the determination of element abundances from the observations as well as in 



the assumed models of chemical evolution of the light element abundances, most importantly 
those of D and 3 He The actual value of the bound, expressed in terms of AN V , has 

therefore become harder to evaluate. Moreover, because the present experimental upper 



bound on the tau-neutrino mass, m Ur < 24 MeV |15|], is relatively close to the upper end of 
the hitherto quoted nucleosynthesis bounds, it would be useful to see how the nucleosynthesis 
can compete with the laboratory when the limit on AN U is considerably weakened. 

The main purposes of this paper are (1) to present a treatment that is accurate enough in 
all the subtleties of computation so that essentially the only uncertainty in the mass bounds 
arises from the abovementioned inherent uncertainty in obtaining neutrino flavor limits from 
matching the SBBN predictions to the observations, and (2) to be general, which is why 
we will present our bounds as functions of the actual nucleosynthesis constraint. We will 
thus always explicitly write down our cross sections and carefully show how we perform 
the thermal averages. Our evolution equations are written in terms of so called pseudo- 



chemical potentials Zi(t) 0, [i7fl , and assume only kinetic equilibrium. This formalism allow 
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us to follow the evolution of the phase space distribution functions instead of the integrated 
number densities and therefore accurately compute the relevant thermodynamic quantities, 
such as the energy - or entropy densities of z/'s. An exception to this is the case of light 
Dirac neutrinos, for which the assumption of kinetic equilibrium does not hold. 

We will include all neutrinos in our equation network. Tracking v e is particularly im- 
portant, because v T v T annihilations below the v e decoupling temperature T Ue ~ 2.3 MeV 
18fl , would produce an excess of u^s around the n/p-freeze out temperature, biasing the 



n *-> p -equilibrium, hence leading to more neutrons being destroyed and therefore to less 
helium being produced. This effect is very large for m u ~ few MeV and affects our bounds 
significantly if AN V > 1. 

In the Dirac case we treat the different helicity populations of the tau neutrino as separate 
species. Again, the physical reason for this is simple: for moderately light tau neutrinos 
the freeze-out temperature is close to the mass scale, so that u^s annihilate while semi- 
relativistic. Due to the chirality of the interaction, positive and negative helicity states 
interact with different strengths at freeze out, and therefore can have different freeze-out 
number densities compared to what one finds when assuming averaged interaction strengths 
and a total equilibrium between helicity populations 0, [5|, [TJ. Somewhat surprisingly, while 
each has a large effect on N u , they compensate each other quite accurately, so as to give final 
total abundance in good accordance with the helicity averaged approach. 

Our final results in the Majorana case agree with some of the earlier results || |7|, when 
restricted to the specific values of the bound on AN U . In the Dirac case we find stronger 
upper limit on the disallowed mass region than ref. @ but weaker than that of ref. ||. 
The upper limits of the excluded regions are particularly sensitive to the changes in AN V , 
opening up a window for a stable tau neutrino below the experimental bound of 24 MeV if 
the nucleosynthesis bound is relaxed to AN v > 1.3 in the Majorana and AN v > 0.8 in the 
Dirac cases respectively. The latter possibility is quite plausible. For the muon neutrino, on 
the other hand, our upper limits can be competitive with the laboratory bound on mass of 
m v < 160 KeV [[0| only for rather restrictive values of AN U < 0.13 in the Majorana, and 
AN U < 0.39 - 0.44 in the Dirac case. 

In section 2 we will derive generic evolution equations for the particle distribution func- 
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tions expressed in terms of the pseudo-chemical potential. We will also discuss some sub- 
tleties of incorporating the time-temperature relationship into the evolution equations. In 
section 3 we discuss how the observational bounds on the helium abundance should be con- 
verted into a bound on N v . In section 4 we consider the Majorana case and in the section 
5 we will derive and solve the equation network with separate equations for the v T helicity 
components in the Dirac case. We will pay special attention to the thermal averaging of the 
helicity amplitudes, which is a nontrivial task because of the lack of the Lorentz invariance of 
the spin dependent matrix elements [pOR . Finally, section 6 contains our conclusions. Some 



calculational details are presented in the appendix. 

2 Generic evolution equations 

In this section we will derive the evolution equations for the particle distribution functions. 
We will also discuss how the time-temperature relation should be consistently incorporated 
to the equation network. Our derivation here relies heavily on that of ref. ]T7[|. We begin by 
writing down a set of Boltzmann equations for the scalar phase space distribution functions: 

E i (d t +pHd p )f t (p,t) = C E4 (p,t) + C^(p,t), (1) 

where E { = (p 2 + m 2 ) 1 / 2 and H = (Mp/ZM^ 1 ' 2 is the Hubble expansion rate, where Mp\ is 
the Planck mass and p is the total energy density. The index % runs over all particle species in 
the plasma; each momentum state in each species has its own equation like ([5]), all of which 
are coupled together through the elastic and inelastic collision terms C#(p, t) and Ci(p,t). 

A tremendous simplification results if the system is in thermal equilibrium; then each 
distribution function can be described by two parameters, the temperature, and possibly, 
a chemical potential. Decoupling particle species however, are by definition not in thermal 
equilibrium. Fortunately though, they often are in close kinetic equilibrium, because the 
kinetic equilibrium is held by elastic scattering processes whose rate typically greatly sur- 
passes that of annihilations, particularly for large m/T. Therefore, at each instant of time, 
the momentum distribution of particles should be closely approximated by a function 

f(p, Zi ) = (e^+* + l)" 1 , (2) 
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where the time dependent function z(t) acts as an effective chemical potential driving the 
system out of the chemical equilibrium. The function z(t) is called pseudo-chemical potential, 
because, unlike the ordinary chemical potential, it appears with the same sign in both the 
distribution function for particles and antiparticles ]TB|. 



We will assume throughout that photons and electrons, because of their extremely fast 
electromagnetic interactions, are in complete thermal equilibrium and therefore we need not 
write down evolution equations for them. For all neutrino species on the other hand it is 
necessary to follow the chemical evolution accurately. In the case of muon and tau neutrinos 
this is obvious, because it is exactly the effect of their energy density on the expansion rate, 
and thereby on the final helium abundance that we wish to study. The electron neutrino is 
known to be nearly massless, so that small changes in the v e number would be of no likely 
importance for the expansion rate. However, even very small variations in n Ve are important 
because of their direct effect weak reaction rates, such as v e + n «-> e~ + p that govern the 
freezeout of these weak interaction rates and the n/p ratio. Using the ansatz (^|) we therefore 
end up with the following equations for the various neutrino pseudo-chemical potentials 



hi + 3Hrii = C a(ij <-> kl) 
{«}< 



(3) 



where i runs through all neutrino species and the sum {a]i is over all the relevant collision 
channels. The compactly written left hand sides of the equations are in fact functions of Zf 

n 

2tt 2 



hi + SHrii = -^lH{Ji{xi,Zi) - xlJ^Xi^Zi)) 



+ ^Ji{xi,Zi) -f y Jo(xi,Zi)^\, (4) 

where the dot refers to time derivative, T 7 is the photon temperature, Tj is the temperature 
of the particle species i, Xi = rrii/Ti, and the functions J n (x, z) are defined as 



e V x 2 +y' 2 +z 



J n (x, z)= dy y\x 2 + y 2 Y 12 y— • (5) 

We cannot write the right hand sides of (^) in terms of the number densities n(zi) unless we 
further approximate the phase space Fermi-Dirac distributions (0) with Maxwell-Boltzmann 



distributions. This additional approximation has often been made in deriving relic abun- 
dances of decoupling species [0, ^ ^1], |22| . We will here keep the more correct FD-statistics 
and postpone writing explicit expressions to the collision integrals to the following sections 
where particular cases are considered. Note however that we dropped the contribution from 
the elastic collision integral, which vanishes under the assumption of kinetic equilibrium. 

The evolution equations (§) are strongly coupled not only through the collision terms, 
but also because of the time-temperature relation; this is particularly explicit in the form (|j) 
for the collision part of the equations (|3|). This complication is a general consequence of the 
assumption of kinetic equilibrium. The usual approach to define the time-temperature rela- 
tionship (which we will find inadequate) is to assume that the energy momentum tensor has 
the particularly simple form T^ v = diag(p, —p — p — p), corresponding to the ideal fluid ap- 
proximation, after which the Einstein equations directly lead into the "energy conservation" 
equation 

p = -3H(p + p), (6) 

where p is the total energy density and p is the total pressure. When energy density and 
pressure are expressed in terms of integrals over particle distributions, equation (||) turns 
into an additional equation relating time and the photon temperature. 

The appearance of other time derivatives, like Tj in (|4]) arises from our choice to param- 
eterize each distribution function by two variables: Zi and Tj. Complete determination of 
the evolution of a system consisting of N separate species would therefore require 2N + 1 
independent equations. It would be possible to obtain additional collision equations to aug- 
ment for example by probing higher moments of the original equation. We will instead 
find it sufficient to follow the simplest physical intuition and assume that the neutrino tem- 
peratures are given by the photon temperature down to the scale where the electrons begin 
to annihilate, and later follow the reference temperature of a completely decoupled massless 
species. That is: 



'4 + 2/i e (T 7; 



1/3 



T Vi ss ^ —^' j T 7 (7) 

where the function h e is related to the electron entropy density by s e = (27r 2 /45)/i e T 7 3 . This 
approach becomes better warranted a posteriori when we find out that the annihilations are 
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always practically complete at temperatures well above the electron annihilation temperature 

^ann — TIl e /3. 

However, even with a well defined closed set of equations, there is a problem with the 



fundamental ideal fluid assumption when dealing with an expanding fluid of a nonrelativistic 
decoupling species. Indeed, it is known [^3[ [24| that in such systems the energy momentum 
tensor acquires new terms such as bulk viscosity. Neglecting these contributions, by sticking 



the energy density in the decoupling species starts to dominate over the rest of the mat- 
ter/radiation in the universe. This only happens at very small temperatures, of course, and 
including bulk viscosity terms would exactly cancel the problematic terms (i^s) in (^). The 
final result of this analysis is that the intuitive approach works well: namely, the photon 
temperature, to a very good approximation evolves as a function of time such that the effect 
of the decoupled species is only felt through their contribution to the total energy density (in 
the Hubble expansion rate). Some straightforward algebra based on this assumption then 
immediately gives the standard formula 



where the function hj is related to the entropy of the interacting species, sj = (27r 2 /45)/i/T 7 3 . 
Combined with equations (0) and (§) the equations @ provides a consistent set of equations 
as the starting point of our analysis. 

3 Nucleosynthesis constraints in terms of AN U 

Let us now outline the procedure that leads to the nucleosynthesis constraints on new particle 
physics models, spelled out in terms of a bound for the effective number of neutrino species 
AN U = N u — 3. The argument goes roughly as follows: whatever the nature of the new 
physical phenomenon, its effect on nucleosynthesis eventually boils down to some calculable 
change in the primordial helium abundance Y4 He . Since the helium abundance on the other 
hand is known to be a monotonic function of energy density of the universe, this change 



direct use of the naive energy conservation law (^). This has to do with the breakdown of the 



to the expression @, eventually leads to a blowup of the time temperature relation when 
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in Ya H( , can be mapped to an effective change in the energy density, which customarily is 
measured in units of energy density corresponding to one massless neutrino species. 

The most stringent nucleosynthesis bounds on arbitrary model parameters are obtained 
if one assumes nothing of the likelihood of the underlying microscopic theory. Consider the 
standard model with N v massless neutrinos as a 'reference theory' which will correspond to 
some unknown extension of standard model. The connection is made at each value of the 
baryon to photon ratio rj = ub/u^, in such a way that in the extended model, the value of 
the 4 He abundance, Y(rj) is matched at the same value of rj to a value of N v in the standard 
model with the same value of Y . This mapping thus has a slight, but eventually negligible 
dependence on rj, for a restricted but relevant range in rj. Next one computes the likelihood 
function for the distribution for the variable N v by comparing BBN predictions with varying 
N v to the data || [10], ||. For example in this was found to lead to the best fit 



N v = 2.2 ±0.3 ±0.4, (9) 

which shows the statistical (from the observational determination of Y and the neutron mean 
life) and systematic uncertainties (from 4 He and to a smaller extent from rj - in (||), it was 
assumed that 7710 = 3.0 ± 0.3). Since one could well imagine theories that would effectively 
lower the value of N v as well as increase it, one has to, in the broadest sense we are discussing 
now, take the bounds @ seriously, and accept that they might show preference for some 
extension of the standard model predicting less helium. Based on this information, the 95% 
CL limit was found to be N v < 3.1 0. 

Systematic errors in the process of inferring the primordial abundances from the obser- 
vations however, are not negligible. The tightest constraints on SBBN for a long time made 
essential use of the inferred upper bound on primordial D+ 3 He-abundance (giving a tight 
lower bound on rj)] this constraint was utilized also in arriving (|9|) ||. It has recently been 
question as to whether or not these abundances are subject to particularly large system- 



atic uncertainties due to their poorly known chemical evolution [njj. Indeed, because both 
chemical and stellar evolution affect the abundances of 3 He, the uncertainty is compounded. 
Standard stellar models predict that low mass stars will be efficient producers of 3 He |[25[| , a 
claim which is seemingly backed up by observations of 3 He in planetary nebulae [EH] . How- 



7 



ever, it appears that when the 3 He yields are included in simple models of galactic chemical 
evolution, no value of 77 leads to concordance with the observed solar and present abundances 
of D and 3 He [p7| . The likely preliminary conclusion is that something is wrong with the 
"standard" models of either chemical and/or stellar evolution as they pertain to 3 He. 

Relying only on the much more robust 4 He and 7 Li abundances leads to a shift downwards 
in the concordance region for rj, and hence to a distribution that peaks much closer to N v = 3. 
Simply taking the observations of 4 He and 7 Li at face value, i.e. without assuming that the 
systematic errors are particularly large to artificially produce concordant values of rj, the 
combined likelihood functions for 4 He and 7 Li show a peak at 7710 = 10 10 77 = 1.8 with a 68% 
CL range of 1.6 - 2.8 and a 95% CL range of 1.4 - 3.8. This range for rj can be translated 



into a most likely value for N v = 2.9. In fact the analog of eq. (w becomes [28 



N v = 3.0 ± 0.3 ± 0.4 (10) 

showing no particular preference to N v < 3 (in fact preferring the standard model result 
of N u = 3) and leading to N v < 4.0 at the 95 % CL level (when adding the errors in 
quadrature). In (|I~0"D, the first set of errors are the statistical uncertainties primarily from 
the observational determination of Y and is identical to the one in (^). The second set of 
errors is the systematic uncertainty arising solely from 4 He, and the last set of errors from 
the uncertainty in the value of 77 and is determined by the combined likelihood functions of 
4 He and 7 Li. 

However, in light of the problems in treating the systematic errors, one might rather take 
a different approach |13 |. Here one assumes the correctness of the standard BBN theory and 



restricts ones scope of extended or modified theories to only those one is deriving bounds 
upon. These new theories have their own prediction of the 4 He abundance, or possibly a 
range of predictions corresponding to the range of acceptable values in their free parameters. 
These new 4 He predictions always correspond to effectively having, say, N v greater than a 
certain critical value N^ mt , the value of which depends on the allowed parameter range in 
the new theory. Thus, for this new theory, all of the distribution in N v below N^ nt must be 
considered unphysical. New, obviously relaxed, bounds on the model parameters follow from 
application of the Bayesian method P3 of cutting the unphysical region of the parameter 
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space and renormalizing the remaining distribution to give approximate (1-a) % CL limits 
on parameters.^ 

In this approach, for the case of stable massive neutrinos, we must use the existing 
laboratory bounds to the extent that i) there are exactly three light neutrinos (LEP) El 
and ii) their masses are further restricted to m v < 160 KeV jB3] and m Ur < 24 MeV []15|| . 
Then in) detailed nucleosynthesis computations show us that within these restricted ranges 
the prediction for helium abundance always corresponds to having N„ > 3, whence the 
unphysical region is determined to be N v < N^ vlt = 3. For example, it has been noted in ref. 
P~3|] that 'strict' bound of N v < 3.13 based on (H), relaxes to a Bayesian bound N v < 3.6 with 



iV£ rit = 3. The more of the distribution lies inside the physical region, the closer the 'strict' 
and Bayesian bounds come to each other. For example the result ( pTOj ) implies a 'strict' bound 
of N u < 4.0 and, to this accuracy, is equivalent to the Bayesian N^ nt = 3 limit. 

After this rather detailed account on how the bounds arise from the nucleosynthesis, we 
wish to stress again that, up to the caveat mentioned in the footnote 2 in case of the Bayesian 
approach, the computation of the nucleosynthesis predictions for a given set of model param- 
eters on one hand, and finding and imposing the observationally derived constraints upon 
them on the other, are unrelated matters. The former can be computed exactly, while one's 
ignorance on the latter can be parameterized with N u . 

4 Majorana case 

We now explicitly develop and solve the evolution equations ([3D for the case in which neu- 
trinos are Majorana particles. We will take the electron neutrino to be massless and let the 
masses of the muon and tau neutrinos vary freely, keeping in mind however, the laboratory 
limits m v < 160 keV ]19[ and vtl Vt < 24 MeV We will assume that electrons and pho- 



tons are in complete thermal equilibrium and write down an equation network comprising 



2 We note, but ignore in the following the slight complication that the bound imposed on model parameters 
by an (1-a) % CL Bayesian limit on N v does not exactly correspond to a Bayesian (1-a) % CL limit imposed 
directly on the parameter space. This is analogous to the case of deducing neutrino mass bounds from decay 
experiments ]2!| , where it is observed that the bound on m is not the same as the root of the bound derived 
for m . 
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all neutrinos. Given the assumptions explained in the previous section, we have 



h Vr + 3Hn Ur = C[y T v T <-> aa) 

h UfX + 3Hn Ufi = Ciy^Vp <-> aa) 

h Ue + 3Hn Ue = Y C[y e v e <-> aa) 



(11) 



where we used the compact notation when writing the left hand side of the equations 
and the dots refer to the equations (|7|) and (§). In practice, we have to isolate the derivatives 



dzj/dT 7 on the left hand side, as (11) is truly a network to solve for the evolution of z^s. It 
would not be practical to show the complicated forms here, however a generic collision term 
appearing on the right hand sides of (|Tl|) is given by: 



J V$ {Zi} I Miypvp -> aa) | 2 S in S fi , (12) 

" spin 

where the symmetry factors Si n and which are equal to unity for the present case, are 
included for completeness. We defined a shorthand notation for the phase space factors 



x 



oo 



dcosfl/ dE p4 1 2 x 
x eWn+^ffa, Zl )f(k 2 , z 2 )f(p 3 , * 3 )/(p4, z 4 ), (13) 

where # is the angle between the incoming 3-momenta ki and k2, is the acoplanarity angle 
between the planes of incoming and outgoing momenta and k = (k\ + k 2 — 2k\k 2 cosO) 1 / 2 . 
The integration limits in the energy of the outgoing particle are 

E^=(E kl+ E k2 ) S + m °- m < TK Xl/2{ \ mlml \ (14) 

max ZS AS 

where X(x, y,z) = (x — y — z) 2 — Axy and s is the usual invariant s = (ki + k 2 ) 2 . The generic 
matrix element is (we always define the momentum labeling as (12 — > 34) in our matrix 



10 



elements) 

J2 \M{vpVp -y aa)\ 2 = Q4G 2 F [(c 2 Va + c\ a ) ((A* • p 4 ) 2 + (k 2 ■ Pa) 2 - ml p 3 ■ Pi) 

spin 

+ (c 2 Va -c 2 Aa )ml (k 1 .k 2 -2ml)}, (15) 

where the normalization of the couplings is such that for neutrinos cy v = ca u = 1/2 and for 
electrons cy e = 2sin 2 w — 1/2 and c^e = —1/2, except in the scattering z/ e z/ e — > ee, where 
due to the additional W-channel, effectively cy e — > 2sin 2 6' w + 1/2 and CA e — > 1/2 after 
a Fierz transformation. While the matrix element (|1^) itself is invariant, the phase space 
distributions are written in the rest frame of the plasma and hence we cannot use the simple 
CM-frame expression for |.M| 2 . Of course, in the approximation where one neglects the final 
state Pauli blocking factors, phase space integral reduces to an 1 dimensional integral over an 
invariant cross section ||22|| . Here we will instead write the dot-products in the frame specified 
in the appendix A and perform the phase space integral without further approximations. 

The numerical solution of (JTTj) proceeds as follows. For a given (pair of) neutrino mass(es) 
we begin with equilibrium distributions at some sufficiently high temperature, in practice at 
Tlnit = 100 MeV, and integrate the equation network down until T en d = 1 KeV (in the photon 
temperature), when nucleosynthesis is essentially over, tabulating the functions Zj(T 7 ) and 
the time temperature relation £(T 7 ) along the integration. Then the resulting interpolation 
tables are used as an input for a properly extended standard nucleosynthesis code, which 
we again run for each mass pair generating isocontours in the primordial helium abundance. 
As described in the previous section, we map the deviation in the helium abundance to a 
deviation in the number of neutrino degrees of freedom: AA^^ = N u — 3. Note that the most 
natural bound is in fact in terms of the helium abundance itself, but we are yielding to what 
has become a the common practice in expressing the nucleosynthesis bounds. 

In figure 1, we plot a specific run showing the change in the electron neutrino density 
due to the annihilation of heavy tau neutrinos with a mass m VT = 5 MeV. Because electron 
neutrinos freeze out at the temperature Tdec(^e) — 2.3 MeV, their number density remains 
close to the equilibrium value until a few MeV. However, since the annihilation of i/ T 's is 
still occurring below that temperature, there is a slight increase in the electron neutrino 
abundance. The excess is about 10% at T 7 = 0.7 MeV, which roughly corresponds to the 

11 



1.25 
1 

f 0.75 
I 0.5 
0.25 


10 1 0.1 0.01 

T 

Figure 1: Shown is the electron neutrino number density corresponding to the v T mass of 
5 MeV and normalized to equilibrium density n = (3C(3)/27r 2 )T 7 3 . The dashed line shows 
for comparison the unperturbed electron neutrino temperature, and the short dashed line 
corresponds to the v T number density. The overall decrease in the densities around T 7 ~ 0.2 
MeV follows from electron-positron annihilations which increase T 7 relative to T Ve . 

temperature when n/p-ratio freezes out. The effect of this excess is to keep n/p-ratio in 
equilibrium until a little later thereby decreasing the amount of neutrons and hence the 
eventual helium abundance; numerically, in the conventional units of AN U this effect is 
found to correspond to ~ —4.65n Ue JHJ, where 8n Ve ~ 0.1 is the actual change in the electron 
neutrino number density (normalized to n eq = 1) in the present example. Combined with the 
opposite effect on the helium abundance due to simultaneous slight increase in the energy 
density, the full effect of the variation in the electron neutrino density in the present example 
is to produce an effective negative contribution of AN^ e ' ~ —0.36 to the number of effective 
species; this example shows the potential importance of accounting for the electron neutrino 
abundance when computing the eventual bounds on masses. 

We have computed the total number of effective neutrino species as a function of neutrino 
mass, either that of or u T , keeping the other neutrinos massless and display the results 
in figure 2. We also show the result for the case where we neglect the effect on the electron 
neutrino density. The effect of electron neutrinos is large in the few MeV region, and it 
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does not affect the eventual bounds for the small values of AN„ very much. However, for 
larger values of AN U the effect can be significant. 

Including the neutrino heating, we find the following bounds on the masses as a function 
of AiV^-bound: in the small mass side 

/MeV < (0.35AiVy 2 + 0.05A^ + 0.59A^/ 2 ) 0(0.15 - AN U ) 

+ (0.09 + 0.47AiV„ + 0.83 ANl - 0.72AJVjJ + 0.26AJV*) 6{AN V - 0.15) 
= <f> M (AN v ), (16) 

where m v is measured in MeVs and 8(x) is the step function. This bound is valid for both the 
muon and the tau neutrino and the error of the fit is less than 1 per cent for 0.01 < AN U < 2. 
In the large mass side: 

mf /MeV > 67.9 - 63. 5 AN v + 38.7AA^ - 15.2AJVj) + 2.4AiV* 

= * M (AN V ), (17) 

which is accurate to better than 1 per cent up to AN U = 2.5. For the both small and large 
mass limits the dependence on 77 is of the order of one per cent for r/ 10 = 1.4 — 3.8. 

For example, using the limit AN U < 1.0 from our equation fllQl) implies the excluded 
region of 0.93 MeV < m* f < 31 MeV. A more stringent bound of AN U < 0.3 would have led 
to bounds 0.31 MeV < m^f < 52 MeV. On the other hand, given the present upper laboratory 



limit of itl Ut < 24 MeV JT^|, opening up a window for a particle in the MeV range would 
require a nucleosynthesis bound weaker than AN V > 1.3. Even with the considerably relaxed 
nucleosynthesis constraints obtained neglecting the D and 3 He data [p8[| , this does not seem 



very likely. Hence the nucleosynthesis bound is still to be viewed very much complementary 
to the laboratory bounds, excluding a stable Majorana neutrino with a mass in excess of 



few hundred KeV. Of course the upper limit found in equation (19) has no relevance for u 



for which the laboratory bound is m v < 160 KeV |BJ. Moreover, the lower bound coming 
form the nucleosynthesis can only be competitive with the laboratory bound if AA^^ < 0.13. 

We complete the Majorana section by noting that the nucleosynthesis bounds are cumu- 
lative; considering the effect of both masses together yields stronger constraints. We have 
computed these bounds by allowing both neutrinos be massive simultaneously in our code. 
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Figure 2: Shown is the effect of a massive neutrino to the nucleosynthesis as a function of 
its mass, expressed in terms of the effective neutrino degrees of freedom. Dotted line for the 
comparison shows the same result ignoring the effect on the electron neutrino density. 

We found that the effect in the <-> v T reaction rates due to their having simultaneously 
nonzero masses is small. Hence the common bounds can be directly derived from (|I1^|T7|). 
For example the low mass limit for m Vr in ([16|) becomes 

m Vr < (f) M {&N v - ^/(m.J). (18) 

A fit for the inverse function (/>m(x), is explicitly given in equation (^) in the section 5.2 
below. Similar expression applies for the large m VT bound with (pM (but not its inverse) 
replaced by $m in (|l8|) . The relative error of these approximate bounds is found to be 

< 5%. 

5 Dirac Case 

In the case of a Dirac neutrino, one is faced with an extra complication resulting from the 
chirality of weak interactions; except in the very nonrelativistic limit, different helicity states 
have different interaction strengths. For m ~ few MeV, neutrinos indeed decouple while 
semi-relativistic, and it behooves us to write independent evolution equations for the two 
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helicity populations. To compare the full treatment to the usual approach which does not 
differentiate between the helicities and using the averaged interaction amplitudes, we note 
the following: First, the R-helicity population interacts much more weakly, and hence their 
relic density gets underestimated in the naive approach. Secondly, the L-helicity population 
interacts more strongly than is assumed in the helicity averaged approach, and there is a 
compensating overestimation of their density. Thirdly, the situation is made more compli- 
cated by existence of t-channel helicity flipping interactions that mix the two species. Clearly, 
in order to obtain high accuracy results, a quantitative computation is required to sort out 
which of these effects is dominant. 

Another problem is that the thermal averaging is more subtle in the Dirac case, because 
the spin dependent matrix elements are not Lorentz invariant. To see this explicitly, consider 
neutrino-neutrino scattering: in naive approach, where one boosts the matrix element to the 
CM-frame, it is (see eg. fl22|)) of the order ~ m^/E 2 . This suppression is particular to the 
CM-frame however, and the true thermal average is in fact of the same order ~ m 2 u as the 
other interactions that dominate in the naive approach [f20 . 



The technical difficulty is greatly increased by the very large number of interaction dia- 
grams, in particular because of a large number of spin flipping t-channel processes that were 
naturally absent in the Majorana case. Finally, the lower bounds on the masses in the Dirac 
case are sensitive to QCD phase transition temperature |2U| because the bounds, as we shall 



see is true even for rather large AN U , are saturated by an out of equilibrium excitation of 
right handed species below Tqcd, the excitation process being the more effective the higher 
Tqcd is. We therefore consider large and small mass cases separately 

5.1 Large mass region 

We first concentrate on the large mass region m > 0(1) MeV. The distinguishing feature 
here is that the particles are heavy enough to have become into equilibrium below the QCD 
phase transition temperature so that their distributions can be described by our kinetic 
equilibrium approach. The region m 0(1) MeV is not well described by the equations 
below and we shall return to this point later (§ |5.2|) . The complete equation network can now 



15 



be written in the following form 



h Ur _ + 3Hn Ur _ = C(v T -U T x <-> aa) + Cg, 



flip 



A=- + 
a=e,v B ,v u 



h Ur+ + 3Hn UT+ = C{v t+ v tX «-> aa) - 



flip 



A=- + 



n Ue + 3Hn Ue = Y C(v e v e <-> aa) 

; , (19) 

where the dots again refer to equations (0) and (§). We included different helicity species 
only for tau neutrinos, because in light of the restrictive laboratory bound on muon neutrino, 
it does not make sense of computing BBN bounds for in the large mass region. The spin 
flip terms appearing ([H]) are given by 



a=e,u e ,u u 



+ { C(v T -V T \ <-> v t +v tX ) + C(z/ T _z/ rA <-> ^r+^rA) } 

A=-,+ 

+ 2C(z/ T _z/ T _ «-> z/ T+ z/ r +), (20) 

where the factor of 2 in the last term accounts for the fact that this interaction changes the 
v— number by 2 units. Generic collision terms appearing in definitions M - are defined 



and normalized similarly to the equations (pT2j- 14) . For example (from now on we will denote 
v T \ by u x ). 

C(u xl u X2 ^aa) = _L^( e *Ai+*A 3 _ e 2*») x 



x 



J V$ {Zi} J^d^Y I M(u xl u X2 -> aa) | 2 , (21) 

^ spin 



where we dropped the symmetry factors which equal to unity and the annihilation matrix 
element is given by 

Y I M(u xl Px2 -> aa) | 2 = 16G^ {(4* + c 2 Aa ) (K X1 ■ p 3 K_ X2 ■ p 4 ) 



spm 



2 
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^(4 a -4J<^Ai-^-A 2 }. (22) 



Here, in order to condense the notation, we dropped the terms that vanish, and combined 
others that are equal under the integration; similar simplifications are made in other matrix 
elements following below. The coefficients cy a and c^a have been defined in section 4 and 
the 4-vector K£ = fc M — m v s^ is related to the 'spin vector' of ith neutrino and can be 
written as 

K£ — (E — \(E 2 - m 2 ) 1/2 )(l; —Xk/k). (23) 

Other matrix elements appearing in the collision terms above include the t-channel scattering 
off the electrons and other neutrinos and their antiparticles. Under the assumption that the 
chemical potentials are small (<^C 1) the distribution functions for particles and antiparticles 
are equivalent, and we can add their contributions under the integral: 

]T | M{v xi f3 -> ^A2/3) | 2 = 16G 2 F Uc 2 Va + c\ a ) (K\i ■ P2 K X1 • p 4 + K X1 ■ Pi K X3 ■ p 2 ) 

spin, /3=a,a 

-(c 2 Va -c 2 Aa )m 2 Q K xl -K X3 }. (24) 
Finally, self scatterings can all be derived from the matrix elements 

I M(is xl v X2 -> u X3 u X4 ) | 2 = 8G 2 F K xl ■ K_ X4 K X2 ■ K_ X3 (25) 

and 

I M(u xl u X2 -> v X3 v XA ) | 2 = 8G 2 F K X1 ■ K X2 K X3 ■ K XA . (26) 

The symmetry factors appearing in the collision integrals are equal to one everywhere except 
the reactions corresponding to (|26|) . There, the symmetry factor is one half, because of the 
degeneracy in either initial or in final state, except for the reaction v T _i> T _ u T+ u T+ , where 
the symmetry factor is 1/4. Note that this reaction changes v T - number by two units, but 
that was explicitly taken into account in the equation (f20|). 

The collision integrals of the massless v e and appearing in ( |T9D are obtained from the 
Majorana matrix element ([15]) in the limit m v = 0. Equations ( |T5| , |22| - F26D exhaust the list of 
interactions relevant for the evolution of the neutrino ensembles. Each of the matrix elements 
( |T5] , P2] - P^D is a polynomial at most of second order of the cosine of the acoplanarity 0. We 
integrate over <fi analytically, after which the remaining 4-dimensional integral is performed 
numerically using the special frame introduced in the appendix A. In figure 3 we show the 
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Figure 3: Shown are the total annihilation rates T Ux± and the total flip rate rg ip (see the 
text) in comparison with the Hubble expansion rate H(T). 

temperature dependence of the annihilation and flip rates T Vx± and r flip , defined by (see 
the appendix A) r, = ^ n j( v M0i cr {ij kl)) for the particular case of m Vr = 5 MeV along 
with the Hubble expansion rate. One sees how the right helicity population drops from the 
equilibrium much earlier than left helicity states. Yet both states are in complete equilibrium 
until well below the QCD phase transition temperature Tqcd ~ few 100 MeV. 

In figure 4, we show a particular example of the evolution of the neutrino energy densities 
as a function of time. While negative and positive helicity populations have different densities 
from each other, their average comes close to that obtained in the helicity averaged approach. 

In figure 5, we show the change in the helium abundance in the Dirac case for >, 1 
MeV. We did our computations also using the helicity averaged approach. The final results 
turned out to be very close to the full solution, in particular in the region of interest for 
nucleosynthesis bounds. While one might have expected this result on qualitatively, the 
quantitative proof only could follow from a numerical calculation. 

We find that the nucleosynthesis bound on the v T mass is fitted to an accuracy of one 
per cent in the range < AN U < 2.5 by: 
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Figure 4: Shown is the evolution of the energy densities of u T _, u T+ and v e corresponding to 
a run of the code with vh Vt = 20 MeV and normalized to po — 77r 2 T 4 /240. Also shown is the 
energy density corresponding to v T in the helicity averaged approximation. 



m°/MeV < 37.8 - 26.9AA^ + 21.3AA^-15.5AA^ + 6.3AA^-0.1AA^ 

= $ D {AN U ). (27) 

This is the main result of this subsection. One observes that the nucleosynthesis constraint 
allows a stable Dirac neutrino just below the laboratory bound m Vr < 24 MeV, given that 
AN V > 0.8. This seems to be admissible given the relaxed constraint following from the 
equation (10|); indeed, the bound AN U < 1.0 gives the constraint > 22 MeV. More 
stringent bound of AiVj, < 0.3 would lead to > 31 MeV. Our result (|2~7l) differs noticeably 
from the previous computations and fall roughly in between the results obtained in || and 
H given the particular values for the bound for AiV,, used therein. 

5.2 Small mass region 

The mass of a neutrino is considered 'small' if the R-helicity population is out of equilibrium 
below Tqcd- Even in this case however, a significant amount of R-helicity states may be 
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Figure 5: Shown is the effective number of degrees of freedom for a Dirac tau neutrino. 



produced by out of equilibrium spin flip processes [ 20f1 . We shall see below that even for very 



large AN U , the lower limit in the exclusion region indeed is saturated by a mass for which 
R-helicity states are out of equilibrium. 

In the small mass region it does not make sense to describe the R-helicity system with 
the distribution function (0). Instead, noting that the R-helicity states are produced from 
an equilibrium ensemble of left helicity states through very mildly energy dependent spin flip 
interactions, the increments in R-helicity population appear with a local (in time) equilibrium 
density characterized by the excitation temperature T(t ex ). Provided that the backscattering 
is not very efficient (to be checked a posteriori) , one can compute the total energy density in 
the right handed species very accurately by including the dilution due to entropy production 



subsequent to excitations. This program was carried out in the reference |20[ . Here we stress 
that the underlying assumption of no backscattering and hence the bounds are valid for 
surprisingly large values of AN V . To this end we extend the treatment of 0] by including 



a simple but accurate model of backscattering. We also point out and correct an inaccurate 
treatment of the effect of a mass of a neutrino for the nucleosynthesis in the final stages of 



the analysis of |2T)fl . The corrected analysis turns out to give bounds roughly 30 per cent 



more stringent than those of ref. p 
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We will write down the starting point of our analysis using the results obtained in ref. 
Because the spin flip processes involving only one right helicity state at the time are by 
far the dominant interactions here, an equation which simply models the back scatterings 
can be written as 

p v + AH Pu = 9ll^LK csT 7 Cu (T) (1 - p/ Peq ) (28) 

where p eq = 77r 2 T 4 /240 is the equilibrium energy density of massless neutrinos, K ~ 16.52 
includes counting over all channels the Fermi-Dirac correction due to the statistics of the 
L-helicity particles and the functions are given by [^(J 

C,(T) = l + 0.81^(T) + 3.71xl0- 2 ^) 2 z 4 ir 2 (z ) 

+(^) 2 4(1 - 2/ 2 ){(l + V 2 )K2{z ± ) - y 2 K (z ± )}) 



C T {T) = 1 + 0.06d r (T) + 3.71 x lO- 2 (^) 2 z *K 2 (z ) (29) 

where the functions di(T) express the nontrivial temperature dependence (deviation from 
the T 7 -law) of the scattering collision term due to the scatterings off muons. The terms 
involving the modified Bessel functions Ki(z) correspond to pion decays with, zq = m^o/T, 
z± = m n ±/T, y = m n ±/m(j, and f w o ~ 93 MeV and f n ± ~ 128.7 MeV. Equation ( |2"8"D is easily 
integrated along with the equation (||) to yield a double integral expression for the relative 
energy density r = p/p eq 

r(x) = /^(^^/^(xOexp f X 'dx"A(x") 

Jx nj\X j Jx 

+ r QCD ( , H ( I{X) , ) 4/3 exp [ X dx'A(x'), (30) 

where tqcd = (17.25/60) 4//3 ~ 0.19 is the diluted energy density of the equilibrium ensemble 
of R-helicity population decoupled above the QCD phase transition, x = T/Tqcd, and 

m = 2.88 ml Qv* (1 + JL%) CM, (31) 

where Tqq D = Tqcd/100 MeV, x = T/Tqq D and m v is in units MeV. Expression (|30| ) obvi- 
ously reduces to those of |2(J when backscattering is neglected. We computed the relevant 
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value of the function r during the nucleosynthesis, r(0), for a large set of parameters m v and 
Tqcd and found that it is to a reasonable approximation 

r(0)^r' QCD e- A ^ + (l-e- A ^), (32) 

where ?"q CD ~ 0.1 and 

Ar M = (2.89 + 2.25T™ D )m^ 

Ar r = (1.29 + 1.34T™ D ) ml (33) 

The accuracy is no worse than 3 per cent for Tqcd = 100 — 200 MeV and AN U < 1.5, which 
corresponds to the whole region of applicability of the final result. For Tqcd = 300 — 400 
it undershoots by about 10 per cent at large AiV^ > 1.5. It should be noted that due to 
the mass effects, a single neutrino with an equilibrium density effectively corresponds to 
1 + f(m) species. Using the nucleosynthesis code we have computed a fit for this function. 
For moderately small masses m u <, 0.6 MeV, it in fact coincides with the function _1 (m y ) 
for the number of effective degrees of freedom for the small mass majorana neutrinos, defined 
in equation (|16|): 

/( m )=0-/( m ) = (-18.6m 3 + 7.9m 2 + 0.02m) 0(0.15 - m) 

+ (0.007m 4 - 0.019m 3 + 0.237m 2 + 1.40m - 0.09) 6(m - 0.15). (34) 

The function f(m) defined above differs considerably from the fit /(m 2 ) used in its place in 
ref. |2O].0 Using the approximation (^) we are finally led to the constraint equation 

f(m) + (1 + /(m))(r QCD e- A ^ + (1 - e~ A ^)) < Atf„. (35) 

We plot the bound for the tau neutrino in figure 6 for Tq C d = 100 — 400 MeV with our 
improved fit function f(m) and using the exact expression ([30]) for r(0). We also show the 
value of r(0) to underline how the effective value for AN U greatly exceeds r(0) for even 
moderate masses. This is the reason why the backscattering correction is so small (we find 
it is typically at most 10 per cent for AN U <, 1.5). 

3 The fit function f(m 2 ) used in p(| unfortunately strongly underestimated the effect of a small neutrino 
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Figure 6: Shown is the bound on the Dirac v T mass as a function of the nucleosynthesis 
constraint AN v for Tqcd — 100 — 400 MeV The two short lines show the function r(0) for 
Tqcd = 100 - 200 MeV. 



There is a slight inaccuracy in the derivation above that goes over the stated approxi- 
mations. Namely, we have computed the mass effect of the excited right handed population 
by multiplying by f(m) their fraction of the energy density r. It is clear however that this 
mass effect actually depends on the relative number density n instead. To correct for this 
inaccuracy would require computing also n in a way similar we found r above. This would 
merely require re-evaluating the collision terms and correcting the power for the dilution fac- 
tor 4/3 — > 1 in equation (^). We will not do so here, because the error made is very small. 
Indeed, using the equation fl32"l) it is easy to show that the relative error (undershoot) in AN U 
is generously bracketed by 6(AN„) < ((a-l)r QCD e- Ar + (6- 1)(1 -e~ Ar ))/(m) < 0.04 AiV„, 
in the mass ranee of interest (a = (60/10. 75) 1/3 and b = (17.25/10.75) 1/3 ). Moreover, this 
error tends to cancel the error made by the approximation (|32|) . 

In case of the muon neutrino, our bound becomes competitive with the laboratory bound 

mass to the nucleosynthesis. This is because it apparently failed to correctly model the dominant source of 
the effect of a neutrino mass for the nucleosynthesis; the change in the capture time of free neutrons. Indeed, 
at the capture temperature T 7 ~ 0.1 MeV, the mass is typically dominating over the radiation, whence one 
expects a strong linear correlation between mass and the induced effective chance in N v , as is seen in our fit 
for f(m) above. 
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of < 160 KeV, when AN U < 0.39 (0.44) for T QCD = 100 (200) MeV. Since we are finding 
that the function f(m) was inaccurately estimated in ref. [|20| we give for comparison the 
correct bounds for the values of Tq CD and AN U < 0.3 discussed there: 

< f 130 KeV T QCD = 100MeV 

~ \ 120 KeV T Q cd = 200 MeV { 1 

< f 150 KeV T QCD = 100MeV 

Ut ~ \ 140 KeV Tqcd = 200 MeV. { 1 

For AN U < 1.0, these limits become: 

< f 310 KeV T Qra = 100MeV 

~ j 290 KeV T QCI) = 200 MeV { 1 

< j 370 KeV T QCD = 100MeV 

^ ~ \ 340 KeV T Q cd = 200 MeV. 1 ; 

Before concluding, we discuss the connection between the computations in the high mass 
and the low mass regions. It should be obvious that when the function r(0) is close to 1, 
one actually enters the region where the kinetic-equilibrium treatment employed at the high 
mass region becomes valid. However, one would expect that the connection of the low and 
high mass solutions is not completely smooth, because close to the crossing point the right 
helicity population is not quite in complete equilibrium, nor completely out of equilibrium 
(tqcd is bigger than assumed in the low mass treatment). Hence it is difficult to improve 
the computation quantitatively without a detailed knowledge of the QCD phase transition 
dynamics. We conclude that the bound fl3"o] ) should be trusted until about AN U ~ 1.0, above 
which there may be large (few tens of per cents) uncontrolled uncertainties in the results. 



6 Conclusions 

In conclusion, we have carefully computed the mass bounds on the stable Majorana and 
Dirac neutrinos arising from nucleosynthesis constraints. We also discussed in detail how 
nucleosynthesis constraints on particle physics models arise, and how (and to what extent) 
they can generically be modeled through the effective number of degrees of freedom. In our 
computation of the mass bounds we included the effects of heating of the electron neutrino 
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system as a result of the annihilations below the v e freeze-out temperature and the effect 
of chirality in the weak interactions on the evolution of different helicity components in 
the case of a Dirac neutrino. We also computed the bounds for the case in which both 
and v T are massive simultaneously, resulting in stronger constraints. Most importantly, we 
computed all bounds as functions of the actual nucleosynthesis constraint on the effective 
number of neutrino degrees of freedom AN U , except in the case of light Dirac neutrinos, 
where nevertheless an implicit function of a form <p D (AN u ,m u , Tq CD ) = was derived for 
the bound. We claim that in all our bounds, the theoretical uncertainty is negligible and 
therefore realistic bounds, or estimates of the uncertainties in the bounds can be obtained 
solely on the basis of a separate analysis of the determination of the bound on AN U . 
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A Phase space integrals 

In this appendix we define the phase space co-ordinate system we employed in our compu- 
tations and write down the matrix element in terms of these variables in a couple examples. 
We took the independent variables to be the magnitude of the two 3-momenta of the in- 
coming states, k\ and k 2 , energy of one of the outgoing states E q2 , the angle between the 
incoming states 9 and the acoplanarity angle <fi between the collision planes. In terms of 
these variables relevant 4-momenta have the expressions (we are using the same expression 
for the 4-momentum and the magnitude of the corresponding 3-momentum; what is meant 
in each occasion should be obvious however) 

ki = (E kl ; k\ sin#i, 0, k\ cos#i) 
k 2 = (E k2 ; -/c 2 sin# 2 , 0, k 2 cos6 2 ) 

P4 = (Ep4; p4 sin 9k cos 0,p4 sin 9k sin 0, p 4 cos9k), (40) 
where the angles are defined as 

cos 6*i = (k\ + k 2 cos 9) / k, 
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cos 9 2 = (k 2 + ki cos9)/k 
1 

kQi 



cos 9 K = — ]— ( E p4 (E kl + E k2 ) - E kl E k2 + hk 2 cos 



— ^(mi + m 2 + m 4 - m 3 ) ) . (41) 



2 

and k = (ki + k 2 ) 2 (cf. equation (|T4|)). In terms of these variables the matrix element ([15]) 
in the Majorana case becomes 

J dcf>\M\ 2 = 1287iG 2 F {(c 2 Va + c 2 Aa )[s 2 /2-2s(E kl E p4 + k 1 p A cos9 1 cos9 K ) 

+AEl 1 E 2 i — SEkiEpJ^iPi cos 9\ cos 9k 
+2k 2 lP 2 4 (2 cos 9\ cos 9 2 K + sin 9\ sin 9 2 K ) - m 2 V(j (s - 2m 2 )) 
+ (c 2 Va -c 2 Aa )m 2 a (s-Gml,)}, (42) 

where s = 2E\ x E k2 — 2k\k 2 cos9 + + m 2 2 . We have checked numerically that when inte- 
grated over the phase space in the Maxwell-Boltzmann approximation, the matrix element 



reproduces the much simpler thermal average over the invariant cross section |22 

1 

512ir 6 nin 2 



(vM0i&) = 7777-7 x f V<£> {0} f d<fr V I M{v p vp -> aa \ 2 S izi S 6 

ol2TT b riTn 2 J Jo 



spm 

MB 1 f°° , r , , 2n ts 



8m 4 TKl{f] 



POO / o 

/ d Sv / ^(s-4m 2 )K 1 (^) ( TcM(s), (43) 

Jim 2 1 



where the cross section is easily obtained by integrating the matrix element ([15] 
<W«) = ^^{(cL + 4j(^(l + ^?)-4m 2 ( S -2m 2 )) 



(44) 



Similar checking is not directly possible in the Dirac case, because there the matrix element 
is not invariant and one does not expect the helicity amplitudes to reduce to the simple 



expression (|43|). However, the helicity summed amplitude is of course again an invariant and 
we undertook to check in every case separately that when summed over initial state helicities 
our thermal averages again did reproduce the simpler results fl4"3"|) over the total scattering 
cross section in the Maxwell-Boltzmann limit. 
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